function diff = transition_t_new( theta, J_end, Q_ini, A, param, aux, ind_t)

    % Note that there are two different Q_t one for the HJB and one for the KFE 
    Q_old=nan( length(aux.t_vec), aux.n_m);
    lnm_grid=nan( length(aux.t_vec), aux.n_m);
    delta_grid=nan( length(aux.t_vec), aux.n_m);
    m_t=repmat([param.m_l_end;param.m_u_end],1,length(aux.t_vec));

    hires_tot   = nan(size(aux.t_vec));
    new_chains  = nan(size(aux.t_vec));
    G_r         = nan(size(aux.t_vec));
    m_l_t       = nan(size(aux.t_vec));
    m_h_t       = nan(size(aux.t_vec));
    m_e_t       = nan(size(aux.t_vec));
    m_u_t       = nan(size(aux.t_vec));

    % Initial guess for boundaries
    if isnan(A(1)) 
        m_l=ones(size(aux.t_vec)).*param.m_l_end;
        m_h=ones(size(aux.t_vec)).*param.m_h_end;
        m_e=ones(size(aux.t_vec)).*param.m_e_end;
    else
        m_l=A(1,:);
        m_h=A(2,:);
        m_e=A(3,:);
    end

    delta_dist  = nan(length(aux.t_vec), length(param.dist_m_grid));
    q_dist      = nan(length(aux.t_vec), length(param.dist_m_grid));


    diff_Q=1;
    nr_Q=1;

    while diff_Q>aux.tol_Q

        excess=nan(size(aux.t_vec));
        mean=nan(size(aux.t_vec));
        mean_m=nan(size(aux.t_vec));
        ete=nan(size(aux.t_vec));

        %% KFE itertation

        % Loop over time from t=0

        for t=1:1:length(aux.t_vec)

            % Calculate lambda on grid
            param.lambda=max(0,interp1(aux.t_lam,theta,aux.t_vec(t),'linear','extrap'));
            param.m_l = m_l(t);
            param.m_h = m_h(t);

            if t==1 && (param.K>0)
                param.m_u=param.m_u_ini*(1-param.shock);
                Q_old(t,:)= Q_ini_fun(Q_ini,param,aux);
            end

            if param.K>0

                param.m_e=m_e(t);

                if t==1
                    tmp_lnm=ln_m_KFE([param.m_l,param.m_u_ini*(1-param.shock)],1,aux);
                    [tmp_delta]=del_num(log(param.m_h), tmp_lnm,Q_old(t,:)',param);            
                else
                    tmp_lnm=ln_m_KFE(m_old,1,aux);
                    [tmp_delta]=del_num(log(param.m_h), tmp_lnm,Q_old(t-1,:)',param);% changed to Q_old
                end


                tmp_d=interp1(tmp_lnm,tmp_delta,log(param.m_e),'pchip'); 

                param.delta_2 = param.m_e.^(-1/(1-param.alpha)).*(tmp_d-(1-param.omega_1)./param.alpha./param.c.*param.m_e+(param.omega_0+param.r.*(param.c+param.K)+param.zeta.*param.c)./param.c);
                tmp_x0= fzero(@ (x) delta_u(abs(x)+param.m_e,param), [0 0.5], optimset('display','off'));

                param.m_u=param.m_e+abs(tmp_x0);

                ln_m=ln_m_KFE([param.m_l,param.m_u],1,aux);

                delta=interp1(tmp_lnm,tmp_delta,ln_m,'pchip');

                delta((ln_m>log(param.m_e)))=delta_u(exp(ln_m((ln_m>log(param.m_e)))),param);
                delta((ln_m>=log(param.m_u)))=0;
                delta(end)=0;
                if max(isnan(delta))
                    error('increase grid')
                end
                clear tmp*
            elseif param.K == 0 
                x0= fzero(@ (x) delta_OJS(abs(x)+param.m_h,param), 10^(-2), optimset('display','off'));
                param.m_u=param.m_h+abs(x0); clear x0

                ln_m=ln_m_KFE([param.m_l,param.m_u],aux)+log(param.m_l);

                delta=param.s.*param.lambda.*ones(size(ln_m));

                delta((exp(ln_m)>=param.m_h))=delta_OJS(exp(ln_m((ln_m>=log(param.m_h)))),param);
                delta((ln_m>=log(param.m_u)))=0;
            end
            m_t(:,t)        = [param.m_l,param.m_u];
            lnm_grid(t,:)   = ln_m;
            delta_grid(t,:) = delta;


            if t==1

                Q_old(t,:)= Q_ini_fun(Q_ini,param,aux);            
                m_old=[param.m_l, param.m_u];

                ln_m_ini = ln_m_KFE(m_old,1,aux);
                q_ini=(Q_old(t,2:end)-Q_old(t,1:end-1))'./(ln_m_ini(2:end)-ln_m_ini(1:end-1));
                mean_ini =((q_ini.*(ln_m_ini(2:end)-ln_m_ini(1:end-1)))'*exp(ln_m_ini(2:end)./(1-param.alpha)));

                n=1;
                aux.dt=aux.t_vec(t)/n;
                aux.n=n;

                ete(t)  = jtj_num( delta, Q_old(t,:)');

                Q_old(t,:)=KFE_t(m_old,Q_old(t,:)',delta,param,aux);         

            else
                aux.n=1;
                aux.dt=(aux.t_vec(t)-aux.t_vec(t-1))./aux.n;
                Q_old(t,:)=KFE_t(m_old,Q_old(t-1,:)',delta,param,aux);
                ete(t)  = jtj_num( delta, Q_old(t-1,:)');
            end    
            mean_m(t) =  mean_t_fun(exp(ln_m),Q_old(t,:)');

            delta_dist( t, :)   = interp1( exp(ln_m), delta, param.dist_m_grid, 'linear','extrap');
            delta_dist( t, :)   = delta_dist(t,:).*(param.dist_m_grid<=exp(ln_m(end))).*(param.dist_m_grid>=exp(ln_m(1))) ...
                                        + param.s.*param.lambda.*(param.dist_m_grid<exp(ln_m(1)));

            q_dist( t, :)   = interp1( exp(ln_m), Q_old(t,:), param.dist_m_grid, 'linear','extrap');
            q_dist( t, :)   = q_dist(t,:).*(param.dist_m_grid<=exp(ln_m(end))).*(param.dist_m_grid>=exp(ln_m(1))) ...
                                + Q_old( t, end).*(param.dist_m_grid>exp(ln_m(end))) ...
                                + Q_old( t, 1).*(param.dist_m_grid<exp(ln_m(1)));

            u_tmp           = Q_old(t,1).*param.s;
            G_r(t)          = interp1( exp(ln_m), (Q_old(t,:)-Q_old(t,1))./(Q_old(t,end)-Q_old(t,1)), param.m_h, 'linear');
            hires_tot(t)    = ete(t).*(1-u_tmp)+param.lambda*u_tmp;
            new_chains(t)   = param.lambda*u_tmp + param.s.*param.lambda.*G_r(t).*(1-u_tmp);

            m_old=[param.m_l, param.m_u];
            q=(Q_old(t,2:end)-Q_old(t,1:end-1))'./(ln_m(2:end)-ln_m(1:end-1));
            mean(t)=((q.*(ln_m(2:end)-ln_m(1:end-1)))'*exp(ln_m(2:end)./(1-param.alpha)));

            %% Calculate a measure of excess demand 
            if t==1
                excess(1)=-(mean(t)-mean_ini) - (mean(t)-mean_ini)/20;
            else
                excess(t)=-(mean(t)-mean(t-1)) - (mean(t)-mean_ini)/20;
            end

        end

        %% HJB iteration

        for t=length(aux.t_vec):-1:1

            if t==1
                aux.dt=aux.t_vec(t);
            else
                aux.dt=aux.t_vec(t)-aux.t_vec(t-1);
            end

            param.lambda=max(0,interp1(aux.t_lam,theta,aux.t_vec(t),'linear','extrap'));
            if t==1
                Q_HJB=interp1(ln_m_KFE(m_t(:,t),1,aux),Q_old(t,:),param.ln_m_HJB,'linear','extrap');        
                Q_HJB(param.ln_m_HJB<log(m_t(1,t)))=Q_old(t,1);
                Q_HJB(param.ln_m_HJB>log(m_t(2,t)))=Q_old(t,end);
            else
                Q_HJB=interp1(ln_m_KFE(m_t(:,t-1),1,aux),Q_old(t-1,:),param.ln_m_HJB,'linear','extrap');        
                Q_HJB(param.ln_m_HJB<log(m_t(1,t-1)))=Q_old(t-1,1);
                Q_HJB(param.ln_m_HJB>log(m_t(2,t-1)))=Q_old(t-1,end);
            end
            if t==length(aux.t_vec)
                J_i=HJB_t(J_end,Q_HJB,param,aux);
            else
                J_i=HJB_t(J_old,Q_HJB,param,aux);
            end
            J_old=J_i;

            m_l_t(t) = interp1( J_old, exp(param.ln_m_HJB), 0, 'pchip');
            m_h_t(t) = interp1( J_old, exp(param.ln_m_HJB), param.c, 'pchip');

            if param.K > 0 
                m_e_t(t) = interp1( J_old, exp(param.ln_m_HJB), param.c + param.K, 'pchip');
                J_old = min(param.c+param.K, max(0,J_old));
            end
        end

        %% Update boundaries
        if nr_Q<5
            aux.adj_q=1;
        elseif nr_Q>10 && nr_Q<=20 
            aux.adj_q=0.1;
        elseif nr_Q>20 
            aux.adj_q=0.01;
        else 
            aux.adj_q=0.5;
        end

        diff_Q = max(max(abs([m_l-m_l_t, m_h-m_h_t, m_e-m_e_t])))
        m_l=m_l+aux.adj_q.*(m_l_t-m_l);
        m_h=m_h+aux.adj_q.*(m_h_t-m_h);
        m_e=m_e+aux.adj_q.*(m_e_t-m_e);
        nr_Q=nr_Q+1;
        %%

    end


    if ind_t == 1

        diff.excess =  interp1(aux.t_vec,excess,aux.t_lam, 'pchip','extrap');
        diff.Q_t=Q_old;

        diff.delta_dist=delta_dist;
        diff.q_dist=q_dist;

        diff.delta=delta_grid;
        diff.m_grid=exp(lnm_grid);
        diff.mean=mean;
        diff.mean_m=mean_m;
        diff.ete = ete;

        diff.lambda=theta;
        diff.t_grid= aux.t_vec;

        diff.G_r=G_r;
        diff.hires_tot= hires_tot;
        diff.new_chains= new_chains;

        diff.m_l_t=m_l_t;
        diff.m_h_t=m_h_t;
        diff.m_e_t=m_e_t;
        diff.m_u_t=m_u_t;

    end    

end